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Abstract. The inverse source problem for the radiative transfer equa- 
tion is considered, with partial data. Here we demonstrate numerical 
computation of the normal operator XyXv where Xv is the partial 
data solution operator to the radiative transfer equation. The numeri- 
cal scheme is based in part on a forward solver designed by F. Monard 
and G. Bal. We will see that one can detect quite well the visible sin- 
gularities of an internal optical source / for generic anisotropic k and a, 
with or without noise added to the accessible data Xv f ■ In particular, 
we use a truncated Neumann series to estimate Xv and X v , which pro- 
vides a good approximation of X v Xy with an error of higher Sobolev 
regularity. This paper provides a visual demonstration of the authors' 
previous work in recovering the microlocally visible singularities of an 
unknown source from partial data. We also give the theoretical analysis 
necessary to justify the smoothness of the remainder when approximat- 
ing the normal operator. 



1. Introduction 

The Radiative Transfer Equation (also sometimes referred to as the linear 
transport or linear Boltzmann equation) is often used to model the propaga- 
tion of particles that exhibit absorption and scattering in various contexts, 
including the behavior of photons within biological tissues or of neutrons in a 
reactor. In this paper, we will consider an inverse problem relevant to optical 
molecular imaging (OMI), which is a physically safe imaging modality that 
has seen some recent advances in research. In this application, biochemical 
markers can be used to detect the presence of specific molecules or genes, 
and suitably designed markers could potentially identify diseases before phe- 
notypical symptoms even appear. Such markers are typically light-emitting 
molecules, such as fluorophores or luminophores, which bind to the desired 
molecule to be detected. In contrast to Single Positron Emission Computed 
Tomography (SPECT), Positron Emission Tomography (PET), or Magnetic 
Resonance Imaging (MRI), optical markers emit low-energy near- infrared 
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photons that are relatively harmless to human tissue. Because of their low- 
energy level, the photons scatter before they are measured. Further specifics 
can be found in the bioengineering literature such as [BJ El [TU EH l32| l2~T]. 

The inverse problem we consider consists of reconstructing the spatial dis- 
tribution of a radiative light source from measurements of outgoing photon 
intensities at the boundary of the medium. In many applications, the prop- 
agation of photons emitted can be modeled as inverse source problems of 
time-independent radiative transfer equations. We will assume that the op- 
tical parameters of the underlying medium are known, so that the problem 
of determining the source is feasible. In practice, such optical parameters do 
not vary too much in biological tissue, and not knowning their true values 
does not significantly effect one's ability to detect edges in the image. It is 
shown in [29j that under mild assumptions on the scattering and absorp- 
tion parameters of the medium this is possible. Other work on the inverse 
source problem for the RTE under varying assumptions can be found in 
[H [261 1221 02] > an d further background on optical tomography can be 
found in [28], Q]. In the case of partial data, [15] shows that one can re- 
cover the visible singularities of the source /, which is a particular subset 
of WF(/). As such, the primary focus of this paper is to provide a numer- 
ical scheme by which to detect such singularities. We now describe more 
precisely the mathematical problem. 

We assume £1 to be a bounded domain in M. n with smooth boundary and 
outer unit normal vector u{x). As in [2l?l [15], we also assume that the data 
is given on the boundary of a larger strictly convex domain Q.\ §> tt. The 
radiative transfer equation on the domain Q with internal source / and zero 
incoming illuminations is given by 



where the absorption a and the collision kernel k are functions with regu- 
larity specified later, the solution u(x, 9) gives the intensity of photons at x 
moving in the direction 9, and d±SQ, is the set of points (x, 9) £ dQ x S n_1 
such that ±v(x) ■ 9 > 0. That is, d±SQ. is the set of points (x, 9) G d£l x S™" 1 
such that 9 is pointing outward or inward, respectively. The source term 
/ will be assumed to depend on x only for our purposes. We also recall, 
given the physical model, that equation 11.11 is only applicable at a single 
frequency, as the parameters a and k typically vary widely with different 
frequencies. It is also important that the photons' energy be relatively low, 
since for high energy photons each scattering event is often accompanied 
by an energy change. This is discussed briefly in [15], and we mention it 
primarily to point out the limitations imposed by the assumptions of the 
model. 





u\a_sn = 
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In the case of full data, we have boundary measurements 
(1-2) Xf(x,6) = u\ a+sn . 

In [29], it is shown that for an open, dense set of absorption and scattering 
coefficients (a, k) G C 2 (Tl x S™" 1 ) x C 2 (Tl x S"" 1 x S™" 1 ), one can recover 
/ G L 2 ($7) uniquely from boundary measurements Xf on all of d+SQ. To 
set up the case of partial data, first let V C d + SQ be open and let V <s V. 
Let xv G Co°(d + SQ,) be a smooth cutoff function such that xv( x ,&) = 1 
for (x, 0) G V and Xv(x, 0) = for (x, 9) ^ V. The boundary measurements 
are then given by 

(1-3) X v f(x, 9) = xv(x, 9)u\ 9+S n- 

To make notation a bit simpler, if V = (complete data) we will just 

write X, since in this case Xy = X. 

We now consider the complementary problem of numerical computation 
of solutions to (II. ip and using the accessible part of the synthethized data 
to visualize the results of |15j . To this end, we employ a technique from 
[19] that uses rotations applied to the parameters in the spectral domain 
to help eliminate the ray effect, which is a byproduct of the discrete ordi- 
nates method. It should be noted that in principle such a method could 
be applied in any dimension, but for the actual computations we will re- 
strict ourselves to the two dimensional case. We also refer the reader to |12] 
for another approach to solving the direct transport equation using finite 
element methods, which we do not use here. 

Our main goal is to approximate the operators X*Xf and XyXyf in 
the two-dimensional case. Ignoring the technical details for the moment, 
the idea is to utilize a Neumann series to approximate the forward solu- 
tion operator Xy and its adjoint Xy. As such, one must truncate the 
series to a finite number of terms when computing XyXy, and the differ- 
ence of the approximation from the true function XyXyf will have some 
(higher) Sobolev regularity. The only theoretical result herein is Lemma 
12.31 which essentially asserts that, if mi terms are computed in the ap- 
proximation of X and mi terms are computed in the approximation of X*, 
then [X*X] approx f - X*Xf G H l+m+1 (n), where m = min{m 1 ,m 2 } and 

feH l (n). 

In Section [2j we review the background related to recovering the visible 
singularities of a source / given partial data of its outgoing intensity on the 
boundary. This is summarized by Theorem 2.2 and Corollary 6.4 of |15j . In 
Section [3] we describe the algorithm used to compute the forward solution 
u to (ll.ip . Following, in Section [4] we detail the similar approach used to 
approximate Xy. Finally, we present some numerical computations for a few 
different types of sources with and without noise added to Xy f in Section 
El Section [5] proves the main technical estimate of this work, and Appendix 
[A] provides some necessary results for iterated weakly singular- type integral 
operators. 
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2. Background and Statement of Main Lemma 
Let us establish some standard notation. We set 

(2.1) T = 9-V x , T 1= T + a, T = T Q + a-K, 

where a denotes the operation of multiplication by o~(x, 9), and K is defined 
by 

(2.2) Kf(x,9)= [ k(x,e,6')f(x,6')d6'. 

In order to compute the normal operator XyXyf, we first need to gen- 
erate the solution to the radiative transport equation (jl.ip . given by 

u = Tf 1 (Id - KT- 1 )- 1 J f 

oo 

= T 1 - 1 Y / (KTr 1 ) m Jf. 

If we set uo(x,0) := [r i ~ 1 /](x, 9) and set ui := T^ x Ku\^\ for / > 1, then 
u = J2i>o u i- I n practice, it is quite simple to compute the scattering term 
Kui^i at each iteration. In order to compute T£~ we will use a simple first 
order Euler method to solve the associated differential equation to which 
T-j - is the solution operator: 

9 ■ V x u(x, 9) + a(x, 9)u(x, 9) = g(x, 9), {x, 9) £ Q, x S™" 1 

u\d-sn = 0. 

For a given source / S L 2 (f2), recall the microlocally visible set corre- 
sponding to partial measurements on d+ST2i, given by 
(2.3) 

M' := {{x,£) G T*n\39 € § n-1 such that 9 • £ = and •)${%, ff) / 0}. 

Here Xy ( x ' ®) ls t ne extension of xv '■ d+SQi -> R to Hi x § n defined by 
X?M) = xy(x + r+(x, 9)9,9). 

We restate the main result of [15] for recovering singularities of a general 
L 2 source with partial data: 

Theorem 2.1 ([15.], Theorem 2.2). Let I be a positive integer. There exists 
an open dense set O t of pairs {a, k) G C°°(ft x S"" 1 ) x C°°(Tl x x S™" 1 x Sg" 1 ) 
depending on I such that given (o~,k) G 0;, «/ (z, £) 6 «W, i/ien there exists 
a function v G H l (Q) such that 

(2.4) (z,0 g W(X^X y /) (z,£) £ W(/ + t;). 



The following corollary applies when H-KT^ || has small enough norm as 
an operator from L 2 (ft x S"" 1 ) -»• L 2 (f} x S™" 1 ). 

Corollary 2.2 ([15j. Corollary 6.4). Suppose that \\KTr \\ < 1. T/ien in 
Theorem \2.1\ we have v = 0. 
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With this in mind, we can state the main theoretical estimate of this work, 
which justifies our approach in the numerical computations. The proof is 
given in Section [5j 

Lemma 2.3. Suppose f G H l (U) and o~,k G C°°. Let mi, ma > 0. Then 

m,2 mi 

j=0 i=0 
where m = min{mi, 777-2} • 

3. Numerical Method for the Direct Problem 

We now present a detailed summary of the method used to solve (jl.ip 
for 77 = 2. Note that at this stage it is not important that the source / be 
isotropic, but such an assumption will be used when computing the normal 
operator later on. As in [19] we use the source iteration method, which 
requires one to solve problems of the form 

9 ■ V x u + a(x, ff)u = g(x, 9) 

(3.1) u\g^ S n = 0- 

Without loss of generality, we may take tt = ID C R 2 where B is the unit 
disk, and assume that a, k and / are all supported compactly in D x S 1 or D 
as appropriate. Such an assumption can be justified by finding a ball B(0, R) 
large enough and rescaling the coordinates accordingly. The main advantage 
here from a numerical standpoint is that Q remains invariant under rota- 
tions. Now, for actually computing Tr 1 , it will be easier to have boundary 
conditions on a cartesian domain. To that end, for each 77 6 (0, 2ir) denote 
9 = 6{rj) = (cos (77), sin (77)) and let 9 = 9(rj) be the counterclockwise ro- 
tation of 9 by 7i/2. That is 9(r/)- L = (— sin 77, cos 77). For each 77 S (0, 27i) 
define the ?7-dependent square 

C v = {x G M 2 such that |x ■ 9\ < 1 and \x ■ ^ < 1}. 

In short, is the square of side length 2 rotated by an angle 77. The 
corresponding incoming and outgoing sets (analogous to d±S£l are given by 

r± jr) = {x G dC v such that x ■ 9 = ±1 and \x ■ 9 L \ < 1}. 

At the heart of the rotational method, we perform a rotational change co- 
ordinates so that the derivative in the direction 9 becomes d x . Fix an angle 
77 and for x = (x,y) G [— 1, l] 2 define 

(3.2) 7x r) (x,y) = [«]^(x,y) := 7i(xcos77 - y sin77,xsin77 + y cos r/, rf) 

(3.3) v v (x, y) = [v] v (x, y) := tj(xcos 77 - y sin 77, xsin 77 + y cos 77) 

where u,v are functions on R 2 and M 2 x S , respectively. The bracket no- 
tation will only be used for functions which already have a subscript. With 
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this change of variables, we can rewrite (|3.1f) as 
d 

(3.4) — u^(x,y) + a v (x,y)u v (x,y) = g v (x,y), 
u v (-l,y) =0 on 

We remark that for a general (possibly non-zero) function h E L 1 (cL SID) , 
we would have the boundary condition h v on obtained by projecting 

h onto T- )V via the relation h r) {x\,X2) = h(PZ (xx,X2,rj)), where 

P± : d±SB B (x,rj) H- P±(x,rj) = ±9(rj) - det(x, O^O^ e r ±)J? . 

Since the measurements are necessarily discrete, we introduce the fixed 
parameters N x , and N scat to represent the number of grid points in each 
spatial dimension, the number of directions measured, and the number of 
scattering terms computed in the series T{~ ^n=o(-^^i~ ) m ^> respectively. 
Typically, we will take each such parameter to be a power of 2, since FFT 
algorithms are used to rotate the grids. Then the basic idea for computing 
T^ 1 is to solve equations of the type f|3.4[) by first computing g v and a v by 
rotating the images of cr(-,rj) and g(-,rj) clockwise by the angle rj in the x 
variable. Then we solve (|3.4p for u v , which can be done, for example, by 
using a simple first order Euler method along each column of the cartesian 
grid. Specifically, denote s x := we set u v (xi,:) = and consider the 
cartesian grid {xi, . . . ,x Nx } <g> {yi, . . .,3^} where x^y, = -1 + (i - \)s x 
for 1 < i < N x . Then the Euler method would give the update 

u v {xj,y) = u v (xj-i,y)+s x {g v {xj-x, y) - cr^(x 3 -_i,y)it t ,(x 3 -_i,y)) , 2 < j < N x . 

In order to compute K, we define the angular step size 8 := and 
sum over the set of angles {ryi, . . . , r/7v d } where r/j = (l — i) 5. We can then 
approximate K by the discrete operator Kg which has the formula 

N d 

(3.5) K s g(x,r]i) = 5^2g(x,r}j)k(x,r}j,r]i), 1 < % < N d . 

3=1 

3.1. Computing Xy . We describe in Algorithm Q] the iterative scheme to 
numerically solve (jl.ip and simulate the partial data Xyf with respect to 
some subset V C 5 + 50. 

What we have not yet described in much detail is the method used to 
compute the rotations of each function, which we discuss briefly within the 
remaining portion of this section. For more details however, we refer the 
reader to [19]. The general idea is to write the rotation map r„(x, y) = 
(x cos rj + y cos rj, — x sin r] + y cos rj) as a composition of dilations and shear- 
ing/slant operations in each variable separately. In particular, we can write 

Trj — 1 ° ^x,— sinr; ° ^y,cos?7 ° ■Sy.tanr/j 

' cos rj 

s y ,a(x,y) = (x,y-ax), s X)/3 (x, y) = (x - /3y, y), q,/3gR 
<Zx,t(x,y) = (tx,y), d y>t (x,y) = (x,ty), i€M. 
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Algorithm 1 Computing Xyf 



Let uq solve 6 ■ Vu + cru = f 
for i = 1 to Nd do 

compute <t„. and f Vi 
solve a ai [t*o] % + a m [u ] Vl = f Vi for [u ] Vi 
«o(x,y,7ft) <- [[«o] w ]_ w (x,y) 
end for 

for j = 1 to iVscat do 
/j «- KgUj-x 
for £ = 1 to iV^ do 
compute 

solve d Xl [uj] m + ^ [ Uj ] Vi = [fj] ni for [nj]^ 
uj(x,y,Vi) <- [[uj] Vi ]_ Vi (x,y) 
end for 
end for 

16: compute xu u 



a+sn 



The shearing/slant operations s Xj/ 3 and s yiQ are implemented in phase space 
using a periodic interpolation function together with some identities utilizing 
the discrete forward and inverse Fourier transform. Specifically, we first 
embed the N x x N x image into a 2N X x N x image, padding the top and 
bottom of the image with ^ x N x arrays of zeros. We then perform the 
vertical shearing operation s y , a , which independently shifts each column of 
the image by an amount that depends linearly on the column index with 
factor —a. 

The operation of shifting a vector x = [x±, ... ,x m ] T by an amount s is 
done in phase space. First we define the 2m-periodic function 

(3.6) D m (y)= } y J yeR. 

msm (i) 

We then define the spectral interpolant 

m 

(3.7) x(y) := ^ xiD m (y - I), 

l=i 

which coincides with Xj when y = j and interpolates between those values 
sinusoidally (see Figure H]). The spectrally shifted vector x s is given by 

x s = [x(l + s),x(2 + s), . . . , x(m + s)] T . 

Note that x s is not a priori defined for s E Z, but it can be extended 
continuously to such points due to the structure of the singularities of D m . 

Now, in practice we have first padded the image above and below with 
zeros, so we have m = 2N X which is even. It is then straightforward to verify 
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Figure 1. The values of a random 16- vector x are plot- 
ted while overlayed with the spectral interpolant x(y) = 
YZiXl D m(y-l)foiye[0,l6}. 



that 

. / .v n N x -1 

g2JV - ft) = 2iVsn^ = 2iV- E ^ ' ^[°' 2 ^ 
2iv x sin 2A ^ /iv a . J= _^ 

Recall the iV-point discrete Fourier and inverse Fourier transforms, given by 

TV 

X(k) = F^ k [x(j)] = Y J x(j)e-^- 1 ^ l \ k = l,...,N, 

1 * 

= Fj^f [X(k)\ = _ £ X{k)e^N{j - l)(k - 1), j = l,...,N. 

k=l 

It is then possible to write x(l — s) for I = 1, . . . ,2^ as a composition 
of discrete Fourier transforms and inverse transforms and multiplication by 
complex exponentials. In particular, 
(3.8) 



The dilation operations d^t and d y j are computed via a resampling done 
in the Fourier domain. In particular, we must resample a vector x of size 
2N X down to a vector x of size N x but with a different step size. We can 
do this by using the iV-point fractional Fourier transform with coefficient a, 
defined by 

JV 

(3.9) X(l) = G^[x} =J2<k)e- 2 ™( k -W-V. 

k=l 

For example, if we start with a vector x sampled at the gridpoints {j — 

and we want to shift x by s (i.e. sample at {j — 1 + and then resample 



-i±(-N x + \){j-l) 
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Angularly Averaged Internal Transport Solution Full Boundary Data Partial Boundary Data 




Figure 2. The angularly averaged solution to (jl.ip corre- 
sponding to the internal source given in Figure [3] (left). The 
other two images show the complete and partial data Xf and 
Xyf, respectively, where V is the outgoing conic set above 
the arc of dD defined by rj £ [0, 7r/3]. The axes are labeled 
such that 9 is the transport direction and rj is the position 
along the boundary <9B, with rj = corresponding to the 
point (1,0). 



back to a vector x taking values at the m points y\ = s + h(l — 1) for 

/ = 1, . . . , m, then 

(3.10) 

S ( y{ ) = J_ e ^(— +i)^~ 1 )G^ [e^ k - m - 1+ ^ s Ff™ k [x^e^^+D^ 
This corresponds to a vertical dilation composed wth a vertical shift (see 

mm)- 



4. Approximating the Normal Operator X v Xy 

After approximating the solution to the forward problem and restrict- 
ing to d+S£l to obtain the simulated data Xf, we can then proceed with 
computing X* Xf, and similarly, XyXyf. Recall that for any N G N 

X v f = xvR+T^Hld - KT^ l )-\jf 

N 

= xvR+T,- 1 (KT^Jf + xvR+T^iKT^f^il - KT^Jf. 
Under the condition that H-fTTj -1 ||^2_ j . i 2 < 1> we then have 

(4.i) x v = j* ( ^([r^fiTTM [xvR+Tr 1 ]* 

\m=0 J 

To numerically compute X v , we separately compute the three types of terms 
appearing in (|4.ip . 



10 MARK HUBENTHAL 

Given that J : L 2 (Q) — > L 2 (Q x S rt_1 ), it is straightforward to verify that 

(4.2) J*g(x) = [ g(x,9)d9. 
We also compute 

(4.3) K*g(x,9)= I k(x,e',9)g(x,9')d6'. 

is™- 1 

In the isotropic scattering case, this gives 

(4.4) K*g(x,9) = k(x) [ g(x,9')d9'. 

We will need the discretized version of (|4.3p in the same vain as ()3.5[) , which 
is given by 

N d 

(4.5) K* & g(x,r]i) = 5^k{x,r] jl r] i )g(x,r] j ) 1 1 < i < N d 

3=1 

One can also verify that Tf 1 : L 2 (Q x S"" 1 ) -> L 2 (Q x S™" 1 ) has adjoint 

(4.6) [T{ l ]*g[x,9) = J exp (- a{x + t9,9) drj g{x + s9,8) ds. 

However, in order to apply [T-f 1 ]* numerically, it is easier to use the asso- 
ciated first order differential equation for which it is the solution operator. 
We already know that T^ 1 is a left inverse for 9 ■ V + a with the boundary 
condition in (jl.ip . So naturally — 9 • V + a is the associated differential 
operator for [T-f 

Thus we can use a first order Euler method to compute [T-," 1 ]* just as with 
[Tf 1 . Finally, we make note of a formula for [xvR+T^ 1 ]* : L 2 (d+S£l, dE) -> 
L 2 (Q, x S n_1 ), which is easily verified to be 

(4.7) [ X vR+Tr l }*g(x, 9) = g*(x, 0) x #(x, 9)E(x, 9). 

Recall the truncation parameter N scat 6 N corresponding to how many 

terms in the Neumann series ^m=o(Pi _1 ]*-^*) m ^° use - We P rocee< i as m 
Algorithm [2] 



5. Smoothness Analysis 

When trying to recover microlocally the most singular part of the source 
via XyXy as analyzed in [15], it turns out that in theory only the first term 
of the Neumann series is needed. This is stated in Lemma 12. 3[ for which 
we now provide the proof, utilizing the weakly singular integral results in 
Appendix lAl 
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Algorithm 2 Computing XyXyf 



v (x,r)) <- [xvR+T^rXyfix^) = (X v f)#(x, V )E(x,r)) 
for j = 1 to iV sca t do 

«- KgVj-i > Apply if* 

for i = 1 to iVrf do 

*%(*iV x ,y) = > Apply [T -1 ]* 

solve -d Xl w Vi + a rh w rh = [vj] Vi for u^. 
^(x,y,?7i) <- w = [w m ]_ v . (x,y) 
end for 
end for 



2 

3: 
4: 
5: 
6: 
7: 
8: 
9: 
10: 
11 



«(x, y) <- 5 Ei=i E;=i at «i(x, y, > Apply J* 



> u is an approximation to XyXyf 



Proof of Lemma \2.SX Observe that 
X* v = [ XV R + T{\1 - KT,- 1 )- 1 J}* 

= Yy Xv R + T^(KT^yj)* + [ XV R + T^{KT^T*+\1 - KT^J)* 

3=0 

X v = xvR+T{\l-KT 1 - l )- 1 J 
mi 

= £ xvR+T^iKT^yj + xvR + T{\KT^) m - +1 (l - KT^J. 
i=0 

Thus 

m,2 mi 

x v x v f - £ ^[ Xy i? + r 1 - 1 (irr 1 - 1 )V]*xvi?+r 1 - 1 (^ 1 - 1 )v/ 

mi 

= J2lXvR+Tr\l - KT^rHKT^r^JYxvR+T^iKT^yjf 

i=0 

m 2 

+ jyxvR + T-\KT^ l y JYxvR+T^iKT^r^Hl - KT^Jf 

3=0 

+ [xvR+TrHl - KT^ 1 y 1 (KT^ 1 ) m2+1 J]*xvR+T^ 1 (KT^ 1 ) mi+1 (I - KT^ l )~ x Jf 
=: Ax +A 2 + A 3 . 
First note that 

[xvR+T^TxvR+T^fix, 0) = xv(x, Of f exp {- a(x + t0, 9) dr^j f(x+s9, 6) ds. 

In particular, [ X v R+T{ 1 ]* xv R+T^ 1 is bounded on x S™" 1 ) for all 

j. Moreover, from the proof of Proposition 6.1 of |15j we have that [(Id — 
KT± can be taken to be bounded on £P(fi x whenever (Id — 

KTi )~ 1 is, which we can assume is the case for j < I + maxjmi, 777,2} + 1. 
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From the discussion in Appendix |Aj we have that 

[(KT r 1 r>+ij]*g(x) = J J Pm2+l{V ^~Zi V ~ *' 0) 9{V, 0) dy d9. 

Using a similar argument as in Lemma 2 of [29] and applying Proposi- 
tion EH we have that [(KT{' 1 ) m2+1 J]* is bounded from H l (Q x S"- 1 ) ->• 

Now, to analyze A\ we write 

mi 

A 1 = [(Krf 1 r 2 +V]*[(Id-/i-Tf 1 )- 1 ]*[ Xl , J R + T 1 - 1 ]*XV'^ + Tf 1 ^(^T 1 - 1 )V, 

i=0 

and from the arguments of the previous paragraph it is clear that A\ maps 
H l (Q) to ^ +m2+1 (0i). Likewise, 

A 2 = YyKT^yjY[ Xv R + T{ l Y Xv R + T{ l {l - KT^ 1 )~ 1 (KT^ l ) mi +1 J 

and so A 2 maps if (fi) to i?' +mi+1 (Oi). A similar argument shows that A3 
maps H\D) to i7'+ mi+m2+2 (^ 1 ). □ 

6. Numerical Computations 

For our numerical computations, the goal is to provide visual verification 
of Theorem 12.11 and more specifically, of the related result of Lemma 12.31 
In all examples, we use a cartesian grid of 256 by 256 with 128 directions 9. 
We'll use the notation (x, y) to denote a point in IR 2 . In order to incorporate 
anisotropy in k, we use the Henyey-Greenstein (H-G) phase function |14j . 
which is very commonly used in optical imaging: 

1-.9 2 „ _ 9 

-■> . ./J Li! 'V ^ J 

v 4^ (l+ 9 2_ 2ff 6i-6»')3/2 ra = 3. 

To set g = corresponds to the isotropic case, while g = 1 and 5 = — 1 
correspond to forward scattering and backscattering, respectively. In typical 
applications g is around 0.80 to 0.95 (|12|). which is characteristic of highly 
forward-peaked scattering. 

We also incorporate noise into the boundary data in the following way. 
The noiseless data Xf is an m x 128 matrix, with the rows corresponding 
to the position along the unit circle and the columns corresponding to the 
particles' direction. Given a parameter [i > 0, for the jth column vj of Xf 

(i.e. [Xf]j) we compute ||«j||2 = J v J v j an d generate a vector w of length m 
with values chosen randomly according to the standard normal distribution 
(variance 1 and mean 0). We then define the jth column of the noisy data 
by 



(6.1) P (e,0') = \ 



(6.2) [Xf noise ]j := [Xf], + [Xf]J [Xf]j —== = v, + ^hj-rr- 

v Vw 1 w \\ w h 
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Source, max = 1 



Normal Operator from Complete Data 



Normal Operator from Partial Data 



• 



I 




Figure 3. A source consisting of circular bump functions 
with height 1. The partial data is measured without noise 
on the set V = V; xamp i es . 



Source, max = 1 



Normal Operator from Complete Data 



Normal Operator from Partial Data 




Figure 4. Same as Figure [3] except with a noise coefficient 
H = 0.50. 

Thus fi represents the fraction of \\vjW2 to which we would like to rescale the 
variance of the noise. 

With these details in mind, for all computed examples we have taken 

I I q gg2 

k(x, y, 9, 9') = — X{x 2 +y 2 <1} (x, y) (0.05 + sin 2 xy) - - - — 2 _ ^ - — - - - 
and 



0"O, V, 0) = 0.5x{ x 2 +y 2 <1} (x, y) [0.05 + cos 2 xy] 



sin 



Moreover, using the notation of Lemma l2.3( we take mi = 8 and m<i = 2. 
This corresponds to computing the first 9 terms of the Neumann series ex- 
pansion for Xf and only 3 terms in the series representation of X* . Finally, 
in all computations we take 

(6.3) V = V; xam pi cs = (r/, 9) G d+SQ = dB x S 1 1 7] e [0, vr/3] and rj ■ 9 > 0}. 



7. Conclusions 

We have presented a numerical method to solve the direct problem for the 
Radiative Transfer Equation (jl.ip based on the work of |19j . which utilizes 
the discrete Fourier transform and fractional discrete Fourier transform to 



Normal Operator from Partial Data 

Figure 5. A source consisting of rectangular bump func- 
tions with height 1. The partial data is measured without 
added noise on the same set V as given in Figure [3l 



Source, max = 1 Normal Operator from Complete Data Normal Operator from Partial Data 




Figure 6. Same as Figure [5] except with a noise coefficient 
/x = 0.5. 

Source, max = 0.99993 Normal Operator from Complete Data Normal Operator from Partial Data 





Figure 7. A spiral pattern of continuous cir- 
cular bump functions of the form g(x, y) = 

Ajl - ^z(x - x ) 2 - jz(y - yo) 2 , where A is the maxi- 
mum height and r is the radius. In this example, starting 
from the largest bump and moving counterclockwise, we 
have the set of heights and radii A = {0.5, 1, 0.3, 0.3, 0.4, 0.3} 
and r = {0.2,0.15,0.1,0.1,0.07,0.03}, respectively. The 
partial data is measured with no added noise on the set V 
given in Figure El 
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Source, max = 0.99993 Normal Operator from Complete Data Normal Operator from Partial Data 




Figure 8. Same as Figure [71 but with an added noise coef- 
ficient fi = 0.50. 

implement a rotation-based method for computing the necessary line inte- 
grals. Moreover, we have computed the normal operator XyXy in a similar 
manner, where Xy is the partial data operator. Ultimately, this has given 
us some nice visual examples in the presence of anisotropic (a, k) with or 
without added noise, where the anisotropic parts of the parameters have 
physically meaningful structure as given in [12} [H] . Such examples visually 
demonstrate the ability to detect the singularities of an unknown source 
with only partial data of the transport solution u at the boundary. 

Most importantly, the results of Figures [3] through [8] show the detection 
of edges resulting from singularities in the respective chosen sources /, and 
this is consistent with what is theoretically predicted in Theorem 12.11 and 
with the microlocally visible set Ai' given in (12. 3ft . Moreover, it provides 
a nice generalization of similar well-known results regarding the detectable 
singularities in the context of limited data computed X-ray tomography, see 
[23]. 
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conversations related to the computational aspects of this work. Much of 
the work in Appendix lAl was completed during a visit to Finland under the 
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Appendix A. Smoothing Properties of Compositions of Weakly 

Singular-type Integral Operators 

Let's consider an operator of type 

where a(x,y,r,r],6) G C c °°(IR n x W 1 x R + x S"" 1 x S™- 1 ). We would like 
to analyze compositions of the form [A m Jf](x, 6) where J : L 2 (R n ) — > 
L 2 (R n x S™ -1 ) is the extension operator Jf(x, 9) = f(x) as used previously. 
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For m = 2, we compute 
[A 2 Jf](x,9) 



-I 
-I 



a(x,y l7 \x - y ± \,x - y 1 ,6) 



ln-1 



[AJf] yi 



x-yi 



x-yi 

a(x,yi, \x - y±\,x - y[,0) y a (z/i)2/2, 



dy 



{x-yil 11 ' 1 J — 3/2 

a(x,yi, \x - y 1 \,x^~y 1 ,6)a (2/1,2/2, \yi - 2/2|, yT^2, z^l/i) 
\x - 2/i| n_1 l2/i - y2\ n ~ 1 



J 



x - yi\ 
yi - 2/2 1 , yi^ r y2, x^y~i 



\n-l 



f(y2)dy 2 dy 1 



dyi 



f(y2)dy 2 - 



Similarly for m = 3, 
[A 3 Jf](x,8) 

a(x,yi, \x - yi\,x^y 1 ,9)a (2/1,2/2, \yi - 2/21, 2/i~^2, z^l/i) 



\x-yi\ n 1 \yi-y2\ n 1 \y2-yz\ n 1 
a (2/2,2/3, 12/2 - 2/31, 2/2^1/3, 2/T-I/2) dy 2 dyi I f(y 3 ) dy 3 



In general, define 
a.\ := a, 

"2(^,2/1,2/2, \x - yi\,\yi - y2\,x^ l ,yT^y 2 ,0) ■= "1(^,2/1, |z - 2/1I, x~^-~yi, 0) 
■ ai (2/1,2/2, 1 2/1 - 2/2|,2/T-^2,^yi) , 

(A.l) 

a m (x,yi, ... ,y m ,\x- 2/1 1 , - • • , \y m -i ~ Um\,x- yi, . . .,y m -i - y m ,0) ■= 

a m - 1(2/1, 2/2, • • • ,2/m, 1 2/1 - 2/2 1 , - • • , \y m -i - ym\,yl - 2/2, • • • ,y m -i -y m ,x - 2/1) 

• a 1 (x,y 1 , \x - yi\,x-y!,0), 

where 2/0 := x. For simplicity of notation, define 

y m := (yi,...,y m -i,y) 
r m ■= (I2/1 - 2/2 1 , - • • , \y m -i - 2/1) 
:= (2/T-I/2, • • • , 2/m^T r 2/)- 



Then set for m > 2 
a m (x,y,9) :-- 



^mi.x , y m , \x yi\,r m ,x y\,v r 
\x - 2/i| n ~ 1 |2/i - 2/2| n_1 • • • \y m -i ~ 



T^jdyi ■■■dym-i- 
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Also let 

T r X , m f a\ a i( x >y> \ x - y\,x-y, 

a := J[5{- - y)\, ai(x,y, 9) := . — = 

\x — y\ n 1 

Then 

[A m Jf](x, 6) = f a m (x, y, 6)f{y) dy, m > 0. 



(A.2) 



Now we would like to show that A m J has a kernel of the form 

P m (x,y, \x - y\,x-y,0) 



\x — y\ n m 



where /3 m (x,y,r,<r),0) is C°°. Ultimately, our goal is to rigorously show for 
m > 2 that A m J is a pseudo differential operator of order — m, by adapting 
the proof of Lemma 2 in which already directly applies to the case 
m = 1. 



A.l. The Case m = 2. For rj € S n_i , define the set V 



l"\(B (l/4)U 



B v (l/4)). Also let A v be the rotational matrix such that A v rj = e±, the first 
unit basis vector of W 1 . Let ip{r) be a smooth, even bump function such that 
< ip < 1, 4> = 1 for \r\ < 1/4, and ip = for \r\ > 1/2. We will frequently 
use the notation r = \x — y\ and r] = x — y. Note that differentiation of 
An with respect to x yields, from the fact that d x x — y = y X , that 

3 x y | 



I , A^ | 



< 



c 



"j '/' — \x— y 

will have bounds of the form 



Similarly, Zth order derivatives of A„ with respect to x 



C 



Of course, if we treat A„ as a matrix 



\x-y\ 1 ' 

of functions of r/ only, and not implicitly dependent on x, then each such 
coordinate function is in C co (S n ~ 1 ). 
We have 



a 2 {x,y, 



a 2 (x,yi,y, \x - yi\, \yi - y\,x - yi,yx - y,0) 



\x - yi\ n 1 - y\ n 1 



dyi. 



We may then cut up the above integral, using ip and suppressing some 
variables, to obtain 



+ 



yi\ 



l - V 



q 2 dyi + 

\x - yi 
\x - y\ 



Iz/i - y\ 
\x - y\ 
y\-y 
\x - y\ 



a 2 dyi 
012 dyi 



=: h + h + h. 

Let's focus on the first term. Substitute w 
y\ = x — \x — y\w and yi — y = \x — y\(x — y — w). We also note that the 



f, so that |ac — y| n <itt; = dyi, 
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region of integration is w € Bq(^). We obtain 



f ^ (^^r) a 2(a;,yi,y, \x - yi\,\yi - y|,^^,yT-l/,0) 

7l = 7 Ix-yil^lyi-yl"- 1 ^ 

k-y|" -2 7 Iw^-^x^y - w] 71 - 1 

(ii- — - ^ x — y — w \ 
x, x — |x — y, |a; — y||u>|, |a; — y||a; — y — w\,w, ,9 dw 
\x-y-w\ I 



r n-2 j L 
= ^=2 J E 



dw. 



ip(\w\)a 2 (x, x - \x — y\w, y, \x — y\\w\, \x — y\\rj — w\,w, rf — w,9) 

^r^n r-=i -dw 

b (i/2) \ w r L \n-w\ n 

ip{\w\)ct2 (x, x — rA~ l w, y,r\w\,r\e\ — w\,w, A~ l e\ — w,9^ 
B (i/2) \w\ n ~ l \e x - w^- 1 

Let's consider the other terms: 

f ^ (^^f) a 2(x,yi,y, \x -yi\,\yi- y\,x^i,y^y,e) 
h = J \x-yi\ n - l \ yi -y\ n - 1 dVl 

4>(\x — y — w\)at2 {x, x — rA~ l w, y, r\w\, r\ei — w\,w, A~ l e\ — w,6^j 



r " 2 7b (i/2) \ w \ n 1 \ e i ~ w \ n 1 



dw 



13 - 7 Ix-j/il^lj/i-yl"- 1 dm 

(l-^(H)-V(lei-H)) 



r n 2 \w\ n 1 1 ei — w\ n 1 



• c*2 [x-i x — rA v 1 w, y, r\w\,r\e\ — w\,w, A v 1 ei — w, 9^j dw 

Notice that after multiplying each term by r n ~ 2 , the remaining integrals are 
smooth in the variables x,y,r,rj and 9. 



A. 2. The General Case. We seek to record a general formula for the inte- 
gral representation of A m Jf{x) which resembles a weakly singular integral 
with integral singularity a ,_^ n _ m . Let's make some new definitions to sim- 
plify the notation. Set yo = x and y m = y m +i = y, and given 1 < j < m — 1, 
define 

w . = x ~yj w > ._ y_i - mi 

3 \x-yj+iY 3 \x-yj+ x \ 
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Note that rj = x — y = w m . Also define wo = for convenience, and then 
set 



Y>(KI) 
Y>(KI) 



We notice a few useful formulas for 1 < j < m — 1 

Wj = wj^i - Wj 

m—l 

yj = x- \x- y\wj ] j \w t \ 
i=j+i 

Vj ~ Vj+i = \x- y j+ i\(w~^i -wj). 

Define 



y := [yi, - ■ ■ ,y m -i,y] 



x — \x 



m—l 

y\wj I 1 \wi\ 

1=3+1 



m—l 



3=1 



w : = 



m—l 

n h 

1=3 



m—l 



3=1 



We can partition the integrations involved in a m (x,y,0) as 



E 



T6{1,2,3}- 



m—l 



a m (x,y, r|w|,u/i 



' ui2— tci ['■'■' \r]— w n 



• i — ^ — — — aw i 

r n-m\ v _ w^n-l J]^ 1 | Wj |n-I-l | Wfc+1 _ 

For < j < m — 1 define 



(iitJi . . . dw 



m—l- 



V), 



B (l/2) 



2* := 5^(1/2) 



Pf := 



«I 3 + lV 

l»\(Bb(l/4)U3ffi— (1/4)). 
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Then for each 7 G {1,2,3}™^, set £> 7 := UT=i V j ■ Then we have the 
integral 



I m (x,y,r,r},6) = ^ / 

7 6{l,2,3}" 1 - 1 ^ 



m— 1 

n ^k-'^^+i - 



(A.3) 



a m (x,y,r|w|,u;i, 



JJ-Wrn-1 



\w2-wi \ ' • • • ' |r/-tu m _i| ' 



r n m ir/-^ m _ii n 1 nnii 1 1^/1" ' 1 nr=i 2 1^+1 - ^i™ 1 

We then define 

(A.4) /3 m (x, y, r, v , 0) := r n ' m I m (x, y, r, rj, 6), 

so that 

P m {x,y, \x - y\,x-y,9) 



dw\ . . . 



A m Jf(x) 



I 



\x - y\ 



-f(y)dy. 



In order to justify the convergence of such an iterated integral (as we 
have changed the order around), we first note that x and y are restricted 
to £1 by assumption. Secondly, the singularities in each variable Wj are 
all locally integrable since near such singularities we have local behavior 
like i — I*-,-, or p=— = 1 — r^— r . It remains to show bounds on the terms in 

\Wj\ n 3 J- \Wj + l-Wj\ n - L 

the sum which involve integrations over one or more of the domains X>| := 
R n \ (S (1/4) U i%-^(l/4)). Let j be the first such index where such 
an integration occurs in a given term. We may integrate out the variables 
w± . . . Wj-i first since the domains of integration are bounded in those cases. 
Since yj € $7, we have x — rwj Yl^Ij+i \ w l\ £ ^- Thus 

diam(O) 



i=i+i m\ 



--: M. 



So we may estimate the integral by 

f C 

n M 

dwj + i • 



n^iiKi"-'- 1 

diam(f]) 2 +^'- n 



< 



u 



n-j-2 



J 1 



r j+2-m 



\W j+2 ~ W j+1 

dw m -i 
C 



n-l 



V, 



Ire— II 



nKiW +1 - 



\V - W m -1\" - W m -2 

- dw j+1 ■ ■ ■ dw m -i, 



n-l 



\W j+2 ~ W j+1 



n-l 



where T>!y is the domain obtained from 2? 7 by removing the domains of 
integration with respect to w± . . . Wj. Since 1 < j < m — 1, we have that 
3 — m < 2 + j — m < 1. In all cases, the behavior in r is an integrable 



NUMERICAL RECOVERY OF SOURCE SINGULARITIES 



21 



singularity since n > 2. Further integrations over domains T>\ will result 
in similar estimates, and in all cases the power on r remains integrable. 
It follows that, upon multiplying the whole integral by r n_1 , we obtain a 
function that is C°° in x, y, r and 9, since the multiplicative factor of r has 
exponent between n — 2 and n + m — 3, all of which are smooth at 0. 

It remains to show smoothness in rj. Notice that when 7 m _i = 1 or 3 (i.e. 
whenu; m _i G5 (l/2) or w m ^ G R ri \(J3o(l/4)U5^(l/4))), then \<r]—w m -i\ 
is uniformly positive. Thus, we can differentiate under the integral sign with 
respect to T] the term 1 ^w=i infinitely many times. On the other hand, 
to differentiate the terms with j m -\ = 2 (i.e. u> m -i G B^l/2)), we make 
the substitution z m _i = rj — w m -\ and note that z m _i G Bq{\/2). Then 

term ] 1 [n _ m = ■, 1 l n ^ m is uniformly bounded away from 0, so we 

may differentiate it arbitrarily many times with respect to r\. The only 
other potential problem will occur if in addition, 7 m _2 = 2 or equivalently, 
w m -2 G B ij j^~ L (l/2) = B i (1/2). In this case, we make yet another 

change of variables z m _ 2 = w^i-w m ~2 = V - z m -\-w m -2, so that z m -i G 
Bq{\/2) and the term ■ — does not depend on ri. Similarly, the 

\W m — l W rn — 2\ 

term , h— , 1 = ■ ■ — is uniformly bounded for z m _o G 

Bq(1/2), and so we can differentiate it arbitrarily many times with respect 
to rj. Continuing in this way, if at any point 7j = 2 we make the change 
of variables Zj = wj+i — Wj and obtain similar bounds. The main idea is 
to transfer the derivative to the numerator function a m . It is a routine 
matter to check that under such change of variables, differentiation of a m 
with respect to rj does not result in any unbounded factors involving rj via 
the chain rule. 

Remark 1. By the same reasoning, if a(x, y, r, 77, 9) is in C^iW 1 xl"x M + x 
x S n_!^ then ^ e cl 

We now consider the following adaptation of Lemma 2 from [TT] : 
Lemma A.l. Let m>l and let A : Co(Q) — > C(r2,S n_1 ) be the operator 

Af(x,9)= [ [ r m ' 1 a{x,r,u,9)f{x + ruo)drduj, 

where a G C°°(ft x R x S"" 1 x S n_1 ). Then A is a classical ^DO of order 
—m with full symbol 

00 

k=m—l 

where 

a k (x,0 = 27r l - [ d^A(x,0,u,9)5 {k) (u-nduj. 
kl J sn -i 
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Proof. The proof directly follows from the proof of the case m = 1 in [11] . 
Let A'(x,r,uj,9) = r m ~ 1 A(x,r,uj,9), and note that if A' is an odd function 
of (r,uj), then Af = 0. So we may replace A 1 = r m ~ 1 A by 

KwM = \{A'{r,u) + A'{-r,-u)) 

= ^(r^Air, u) + (-r) m - x A(-r, -uj)) 
r m—l 

= —-{A{r^) + (-ir- 1 A(-r, -uj)). 
We can then integrate over r > only and double the result. Thus, 

P f'OO 

Af(x) = 2 / A' cvcn (x,r,uj,9)f(x + ruj)drduj. 

Js™- 1 Jo 

Changing to polar coordinates via z = ruo and setting y = x + rw, we obtain 

Af{x) = 2 J A' cvcn (xAy-xlV— E- o\ [ ^ dy. 

We then use a finite Taylor expansion of A even (x, r, uj, 9) in r near r = 
with N > to get 

N-1 

A' even {x,r,uj,9) = ^^C^fcCx.w, 0)r k + r N R N (x,r,u,9). 

k=0 

One can compute that 2vlg ven fc (x, w, #) = A' fc (ic, w, 6*) + (— l) fc ^(x, — uj, 9), 
where = d k \ r=0 A' = d k \ r L [r m ~ l A] . Clearly, A' k = for all k < m-1. 
Therefore, following the proof of Lemma 2 in |11| . we obtain that the terms 
a k ( x , £, 9) = 2ni k A' k (x, uj, 9)5^ (u ■ £) dw = for all k < m - 1. □ 

Proposition A.l. Let a £ C£°(R n X R n X R X S"" 1 x S"" 1 ) and consider 
the operator A : L 2 (R™ x S"" 1 ) -> L 2 (R"; C^S"" 1 )) de/iraed 6y 

m<*.»> - / ° (x,l,, |L': y 'iLv r "'' ,6l) /(».»' r »)*- 

Then for m > 1, ^4 m J is a classical pseudodifferential operator of order 
—m with smooth parameter 9, where J : L 2 (M n ) — > L 2 (R™ x S n_1 ) is the 
extension operator Jf{x,9) = f(x). 
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Proof. Recall that 

[A m Jf](x, 9) = I I m (x, y, \x — y\,x -y, 9)j[j,) ,],, 
Pm(x,y, \x - y\,x~^y, 



\x — y\ n m 



-f(y)dy 



I 

p poo 

/ / r m_1 j3 m (x,x — ruj,r,oj,9)f(x — ruo) dr du 
Jgn-i Jo 

r m ~ 1 f3 m (x, x + ru, r, —uj, 9)f(x + ru>) dr doj. 



We can then apply Lemma lA.ll to complete proof. □ 

Let us also compute the adjoint of the operator A m J. Given functions 
/ G L 2 (R n x S™- 1 ) and g G L 2 (R n ), we have 

([A m J}*f,g) LHRn) 

= (f,A m Jg) L 2( R n x§ n-l* ) 

Pm(x,y, \x - y\,x^~y,e) 




x — y\ n m 



-9(y)f(x,0)dydxd9 



Mx,y,\x-y\,x-y,9)_ f ^ d)dxd \ ^ 



\x — y \ n m 



Thus 



[A m jrf(x)=[ [ ^^^ r ^f(y,9)d9dy. 

J R nJ S n-l \X-y\ n m 

It is then clear that [A m J]* : H l (n x S"- 1 ) -)■ H l+m {Vt) using a similar 
argument as in Lemma 2 of [29] together with Proposition I A . 1 [ 
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